#include <src/pcps.h>
#include <src/input.h>
#include <src/operator.h>
#include <src/psi3_interface.h>
#include <src/molpro_interface.h>
#include <src/hamiltonian.h>

///////////////////////////////////////////////////////////////////////////////////////////////////
/// \brief   creates text files with the one and two electron integrals in dmrg input format
///
/// \param[in]   userinp  user input options
///
///////////////////////////////////////////////////////////////////////////////////////////////////

template <class SCALAR> void pcps::Hamiltonian<SCALAR>::write_dmrg_ints(const pcps::Input &userinp) const {

  // get the file names for the one and two electron integrals
  const std::string int1_file = userinp.workdir() + "dmrg_int1.txt";
  const std::string int2_file = userinp.workdir() + "dmrg_int2.txt";

  // open the one electron integral file
  std::ofstream ofs_int1(int1_file.c_str());

  // check that the file is open
  if (!ofs_int1.is_open())
    throw pcps::Exception("failed to open file \"%s\" in pcps::Hamiltonian::write_dmrg_ints") % int1_file;

  // write the number of orbitals
  ofs_int1 << _norb << std::endl;

  // write the one electron integrals
  for (int p = 0; p < _norb; p++)
  for (int q = 0; q < _norb; q++)
    ofs_int1 << boost::format("%4i  %4i  %60s\n") % p % q % pcps::format_number("%22.12e", this->get_oei(p, q));

  // close the one electron integral file
  ofs_int1.close();

  // open the two electron integral file
  std::ofstream ofs_int2(int2_file.c_str());

  // check that the file is open
  if (!ofs_int2.is_open())
    throw pcps::Exception("failed to open file \"%s\" in pcps::Hamiltonian::write_dmrg_ints") % int2_file;

  // write the number of orbitals
  ofs_int2 << _norb << std::endl;

  // write the two electron integrals
  for (int p = 0; p < _norb; p++)
  for (int q = 0; q < _norb; q++)
  for (int r = 0; r < _norb; r++)
  for (int s = 0; s < _norb; s++)
    ofs_int2 << boost::format("%4i  %4i  %4i  %4i  %60s\n") % p % r % q % s % pcps::format_number("%22.12e", this->get_tei(p, q, r, s));

  // close the two electron integral file
  ofs_int2.close();

}

// explicitly instantiate the function
template void pcps::Hamiltonian< double               >::write_dmrg_ints(const pcps::Input &) const;
template void pcps::Hamiltonian< std::complex<double> >::write_dmrg_ints(const pcps::Input &) const;

///////////////////////////////////////////////////////////////////////////////////////////////////
/// \brief   prepares the Hubbard model Hamiltonian for fermions or spinless fermions
///
/// \param[in]   userinp  user input options
///
///////////////////////////////////////////////////////////////////////////////////////////////////

template <class SCALAR> void pcps::Hamiltonian<SCALAR>::initialize_hubbard(const pcps::Input &userinp) {

  // get the model parameters
  const SCALAR t = pcps::unity<SCALAR>() * userinp.hubbard_t();
  const SCALAR u = pcps::unity<SCALAR>() * userinp.hubbard_u();

  // get all pairs of nearest neighbor sites
  const std::set< std::pair<int,int> > nn_pairs = pcps::get_nn_site_pairs(userinp);

  // spinless fermion hubbard model
  if (userinp.particle_type() == pcps::keyword::Part_Spinless_Fermion ) {

    // create operators for each nearest neighbor pair
    for (std::set< std::pair<int,int> >::const_iterator pair = nn_pairs.begin(); pair != nn_pairs.end(); pair++) {

      // create hopping operators
      { OpBaseConstPtr temp_ptr( new pcps::OpRL<SCALAR>(pair->first, pair->second, -t) ); _ops.push_back(temp_ptr); }
      { OpBaseConstPtr temp_ptr( new pcps::OpRL<SCALAR>(pair->second, pair->first, -t) ); _ops.push_back(temp_ptr); }

      // create repulsion operators
      { OpBaseConstPtr temp_ptr( new pcps::OpNN<SCALAR>(pair->first, pair->second,  u) ); _ops.push_back(temp_ptr); }

    }

  // full fermion hubbard model
  } else if (userinp.particle_type() == pcps::keyword::Part_Fermion ) {

    // create operators for each nearest neighbor pair
    for (std::set< std::pair<int,int> >::const_iterator pair = nn_pairs.begin(); pair != nn_pairs.end(); pair++) {

      // create alpha hopping operators
      { OpBaseConstPtr temp_ptr( new pcps::OpRL<SCALAR>(2*pair->first, 2*pair->second, -t) ); _ops.push_back(temp_ptr); }
      { OpBaseConstPtr temp_ptr( new pcps::OpRL<SCALAR>(2*pair->second, 2*pair->first, -t) ); _ops.push_back(temp_ptr); }

      // create beta hopping operators
      { OpBaseConstPtr temp_ptr( new pcps::OpRL<SCALAR>(2*pair->first+1, 2*pair->second+1, -t) ); _ops.push_back(temp_ptr); }
      { OpBaseConstPtr temp_ptr( new pcps::OpRL<SCALAR>(2*pair->second+1, 2*pair->first+1, -t) ); _ops.push_back(temp_ptr); }

    }

    // create onsite repulsion operators
    for (int i = 0; i < userinp.nsites() / 2; i++)
      { OpBaseConstPtr temp_ptr( new pcps::OpNN<SCALAR>(2*i, 2*i+1, u) ); _ops.push_back(temp_ptr); }

  // unknown particle type
  } else
    throw pcps::Exception( (boost::format("Hubbard model not defined for particles of type \"%s\"") % userinp.particle_type()).str() );

}

///////////////////////////////////////////////////////////////////////////////////////////////////
/// \brief   prepares the Heisenberg model Hamiltonian
///
/// \param[in]   userinp  user input options
///
///////////////////////////////////////////////////////////////////////////////////////////////////

template <class SCALAR> void pcps::Hamiltonian<SCALAR>::initialize_heisenberg(const pcps::Input &userinp) {

  // check that the particles are spins
  if (userinp.particle_type() != pcps::keyword::Part_Spin)
    throw pcps::Exception( (boost::format("heisenberg model requires the particles to be of type \"%s\"") % pcps::keyword::Part_Spin).str() );

  // get the particles' spin
  const double s = userinp.particle_spin();

  // get the external field strength (units of J)
  const SCALAR b = pcps::unity<SCALAR>() * userinp.heisenberg_b();

  // create vectors to hold diagonal and off-diagonal operators
  std::vector<OpBaseConstPtr> diagonal_ops;
  std::vector<OpBaseConstPtr> offdiagl_ops;

  // get a set of all nearest neighbor pairs
  const std::set< std::pair<int,int> > nn_pairs = pcps::get_nn_site_pairs(userinp);

  // initialize a vector to hold the coordination number of each site
  std::vector<int> coords(userinp.nsites(), 0);

  // process each nearest neighbor pair
  for (std::set< std::pair<int,int> >::const_iterator pair = nn_pairs.begin(); pair != nn_pairs.end(); pair++) {

    // create short names for the two sites
    const int i = pair->first;
    const int j = pair->second;
    assert(i != j);
    assert(i < userinp.nsites());
    assert(j < userinp.nsites());

    //// add imaginary time operator (assumes spin-1/2 spins)
    //if (userinp.optimization() == pcps::keyword::Opt_Imaginary_Time) {
    //  if (std::abs(userinp.particle_spin() - 0.5) > 1.0e-9)
    //    throw pcps::Exception("imaginary time evolution only implemented for spin-1/2 spins");
    //  if (std::abs(b) > 1.0e-9)
    //    throw pcps::Exception("imaginary time evolution only implemented for zero magnetic field");
    //  _imag_time_sites.push_back(std::vector<int>(2));
    //  _imag_time_sites.rbegin()->at(0) = i;
    //  _imag_time_sites.rbegin()->at(1) = j;
    //  _imag_time_ops.push_back(std::vector<double>(4*4, 0.0));  //  down down = 0, down up = 1, up down = 2, up up = 3
    //  _imag_time_ops.rbegin()->at(0*4+0) =  0.25;
    //  _imag_time_ops.rbegin()->at(1*4+1) = -0.25;
    //  _imag_time_ops.rbegin()->at(2*4+2) = -0.25;
    //  _imag_time_ops.rbegin()->at(3*4+3) =  0.25;
    //  _imag_time_ops.rbegin()->at(1*4+2) =  0.50;
    //  _imag_time_ops.rbegin()->at(2*4+1) =  0.50;
    //}

    // add raising/lowering operators
    { OpBaseConstPtr temp_ptr(new pcps::OpRL<SCALAR>(i,j,0.5*pcps::unity<SCALAR>())); offdiagl_ops.push_back(temp_ptr); }
    { OpBaseConstPtr temp_ptr(new pcps::OpRL<SCALAR>(j,i,0.5*pcps::unity<SCALAR>())); offdiagl_ops.push_back(temp_ptr); }

    // add S_z operator
    { OpBaseConstPtr temp_ptr(new pcps::OpNN<SCALAR>(i,j,1.0*pcps::unity<SCALAR>())); diagonal_ops.push_back(temp_ptr); }

    // count the sites' coordination
    coords.at(i)++;
    coords.at(j)++;

  }

  // add single-site number operators
  for (int i = 0; i < coords.size(); i++)
    if (coords.at(i) > 0 || std::abs(b) > 1.0e-9)
      { OpBaseConstPtr temp_ptr( new pcps::OpN<SCALAR>(i, b - pcps::unity<SCALAR>() * s * double(coords.at(i))) ); diagonal_ops.push_back(temp_ptr); }

  // add constant offset operator
  if (nn_pairs.size() > 0 || std::abs(b) > 1.0e-9) {
    OpBaseConstPtr temp_ptr( new pcps::OpIdentity<SCALAR>(pcps::unity<SCALAR>() * s * s * double(nn_pairs.size()) - userinp.nsites() * s * b) );
    diagonal_ops.push_back(temp_ptr);
  }

  // place the operators in the hamiltonian so as to optimize parallelism
  const double ratio = double( offdiagl_ops.size() ) / double ( diagonal_ops.size() );
  int j = 0;
  for (int i = 0; i < diagonal_ops.size(); i++) {
    _ops.push_back(diagonal_ops.at(i));
    for ( ; j < ratio * i; j++)
      _ops.push_back(offdiagl_ops.at(j));
  }
  for ( ; j < offdiagl_ops.size(); j++)
    _ops.push_back(offdiagl_ops.at(j));

}

///////////////////////////////////////////////////////////////////////////////////////////////////
/// \brief   returns a set containing all pairs of sites
///
/// \param[in]   nsites  the number of lattice sites
///
/// \return  a set of all possible pairs of sites
///
///////////////////////////////////////////////////////////////////////////////////////////////////

std::set< std::pair<int,int> > pcps::get_all_site_pairs(const int nsites) {

  // initialize return value
  std::set< std::pair<int,int> > retval;

  // add each unique site pair to the set
  for (int i = 0; i < nsites; i++)
    for (int j = i+1; j < nsites; j++)
      retval.insert(std::pair<int,int>(i,j));

  // return the set of all site pairs
  return retval;

}

///////////////////////////////////////////////////////////////////////////////////////////////////
/// \brief   returns a set containing all nearest neighbor site pairs
///
/// \param[in]   userinp  user input options
///
/// \return  a set of all pairs of nearest neighbor sites
///
///////////////////////////////////////////////////////////////////////////////////////////////////

std::set< std::pair<int,int> > pcps::get_nn_site_pairs(const pcps::Input &userinp) {

  // initialize the return value
  std::set< std::pair<int,int> > retval;

  // cartesian and triangular lattices
  if (userinp.lattice_type() == pcps::keyword::Lattice_Cartesian || userinp.lattice_type() == pcps::keyword::Lattice_Triangular) {

    // get the lattice parameters
    const int length_x = userinp.lattice_length_x();
    const int length_y = userinp.lattice_length_y();
    const int length_z = userinp.lattice_length_z();
    const bool periodic = ( userinp.boundary() == pcps::keyword::Boundary_Periodic );
    const bool triangular = ( userinp.lattice_type() == pcps::keyword::Lattice_Triangular );

    // check that a triangular lattice is 2 dimensional
    if (triangular && length_z != 1)
      throw pcps::Exception("lattice_length_z must equal 1 for the triangular lattice");

    // check that no lengths are negative
    if (length_x <= 0 || length_y <= 0 || length_z <= 0)
      throw pcps::Exception("pcps::get_nn_site_pairs requires lattice lengths to be positive integers");

    // loop over all sites in the lattice
    for (int i = 0; i < length_x; i++) 
    for (int j = 0; j < length_y; j++) 
    for (int k = 0; k < length_z; k++) {

      // if the site has a nearest neighbor one position above it along the x axis, add that pair
      if (periodic || i < length_x - 1) {
        const int I = length_x * length_y * k + length_x * j + i;
        const int J = length_x * length_y * k + length_x * j + (i+1) % length_x;
        if (I != J)
          retval.insert( std::pair<int,int>(std::min(I,J), std::max(I,J)) );
      }

      // if the site has a nearest neighbor one position above it along the y axis, add that pair
      if (periodic || j < length_y - 1) {
        const int I = length_x * length_y * k + length_x * j                    + i;
        const int J = length_x * length_y * k + length_x * ( (j+1) % length_y ) + i;
        if (I != J)
          retval.insert( std::pair<int,int>(std::min(I,J), std::max(I,J)) );
      }

      // if the site has a nearest neighbor one position above it along the z axis, add that pair
      if (periodic || k < length_z - 1) {
        const int I = length_x * length_y * k                    + length_x * j + i;
        const int J = length_x * length_y * ( (k+1) % length_z ) + length_x * j + i;
        if (I != J)
          retval.insert( std::pair<int,int>(std::min(I,J), std::max(I,J)) );
      }

      // if the lattice is triangular and has a nearest neighbor up and to the right, add that pair
      if ( triangular && ( periodic || (i < length_x - 1 && j < length_y - 1) ) ) {
        const int I = length_x * length_y * k + length_x * j + i;
        const int J = length_x * length_y * k + length_x * ( (j+1) % length_y ) + (i+1) % length_x;
        if (I != J)
          retval.insert( std::pair<int,int>(std::min(I,J), std::max(I,J)) );
      }

    }

  // icosidodecahedron lattice
  } else if (userinp.lattice_type() == pcps::keyword::Lattice_Icosidodecahedron) {

    // get a set of all 20 triangles on the lattice
    std::set< std::set<int> > triangles;
    {
      std::set<int> s;
      s.clear(); s.insert(24-1); s.insert(27-1); s.insert(30-1); triangles.insert(s);
      s.clear(); s.insert(30-1); s.insert( 3-1); s.insert( 6-1); triangles.insert(s);
      s.clear(); s.insert( 6-1); s.insert( 9-1); s.insert(12-1); triangles.insert(s);
      s.clear(); s.insert(12-1); s.insert(15-1); s.insert(18-1); triangles.insert(s);
      s.clear(); s.insert(18-1); s.insert(21-1); s.insert(24-1); triangles.insert(s);
      s.clear(); s.insert( 1-1); s.insert( 3-1); s.insert( 4-1); triangles.insert(s);
      s.clear(); s.insert( 7-1); s.insert( 9-1); s.insert(10-1); triangles.insert(s);
      s.clear(); s.insert(13-1); s.insert(15-1); s.insert(16-1); triangles.insert(s);
      s.clear(); s.insert(19-1); s.insert(21-1); s.insert(22-1); triangles.insert(s);
      s.clear(); s.insert(25-1); s.insert(27-1); s.insert(28-1); triangles.insert(s);
      s.clear(); s.insert( 4-1); s.insert( 5-1); s.insert( 7-1); triangles.insert(s);
      s.clear(); s.insert(10-1); s.insert(11-1); s.insert(13-1); triangles.insert(s);
      s.clear(); s.insert(16-1); s.insert(17-1); s.insert(19-1); triangles.insert(s);
      s.clear(); s.insert(22-1); s.insert(23-1); s.insert(25-1); triangles.insert(s);
      s.clear(); s.insert(28-1); s.insert(29-1); s.insert( 1-1); triangles.insert(s);
      s.clear(); s.insert( 2-1); s.insert( 5-1); s.insert( 8-1); triangles.insert(s);
      s.clear(); s.insert( 8-1); s.insert(11-1); s.insert(14-1); triangles.insert(s);
      s.clear(); s.insert(14-1); s.insert(17-1); s.insert(20-1); triangles.insert(s);
      s.clear(); s.insert(20-1); s.insert(23-1); s.insert(26-1); triangles.insert(s);
      s.clear(); s.insert(26-1); s.insert(29-1); s.insert( 2-1); triangles.insert(s);
    }
    assert(triangles.size() == 20);

    // for each triangle, create its three nearest neighbor pairs
    for (std::set< std::set<int> >::const_iterator t = triangles.begin(); t != triangles.end(); t++) {
      assert(t->size() == 3);
      std::set<int>::const_iterator i = t->begin();
      std::set<int>::const_iterator j = t->begin();
      j++; assert(*i >= 0 && *j > *i); retval.insert( std::pair<int,int>(*i, *j) );
      j++; assert(*i >= 0 && *j > *i); retval.insert( std::pair<int,int>(*i, *j) );
      i++; assert(*i >= 0 && *j > *i); retval.insert( std::pair<int,int>(*i, *j) );
    }
    assert(retval.size() == 60);

  // unknown lattice type
  } else
    throw pcps::Exception( (boost::format("error in pcps::get_nn_site_pairs:  unknown lattice type \"%s\"") % userinp.lattice_type()).str() );

  // check that the pairs are ordered and do not contain duplicate indices
  for (std::set< std::pair<int,int> >::const_iterator it = retval.begin(); it != retval.end(); it++) {
    assert(it->first < it->second);
  }

  // return the set of all nearest neighbor pairs
  return retval;

}

///////////////////////////////////////////////////////////////////////////////////////////////////
/// \brief   converts the words of an input file line into a hamiltonian operator
///
/// \param[in]   line  a line from the Hamiltonian text file
///
///////////////////////////////////////////////////////////////////////////////////////////////////

template <class SCALAR> void pcps::Hamiltonian<SCALAR>::add_operator_from_file(const std::string & line) {

  // read the line using a string stream
  std::stringstream ss(line);

  // get the keyword
  std::string key;
  ss >> key;

  SCALAR x;
  int i, j, k, l;

  // do nothing if the line is blank
  if (key.size() == 0) return;

  // add an identity operator
  else if (key == "i") {
    ss >> x;
    OpBaseConstPtr temp_ptr( new pcps::OpIdentity<SCALAR>(x) );
    _ops.push_back(temp_ptr);
  }

  // add a number operator
  else if (key == "n") {
    ss >> i >> x;
    OpBaseConstPtr temp_ptr( new pcps::OpN<SCALAR>(i, x) );
    _ops.push_back(temp_ptr);
  }

  // add a double number operator
  else if (key == "nn") {
    ss >> i >> j >> x;
    OpBaseConstPtr temp_ptr( new pcps::OpNN<SCALAR>(i, j, x) );
    _ops.push_back(temp_ptr);
  }

  // add a raising-lowering operator
  else if (key == "rl") {
    ss >> i >> j >> x;
    OpBaseConstPtr temp_ptr( new pcps::OpRL<SCALAR>(i, j, x) );
    _ops.push_back(temp_ptr);
  }

  // add a double raising-lowering operator
  else if (key == "rlrl") {
    ss >> i >> j >> k >> l >> x;
    OpBaseConstPtr temp_ptr( new pcps::OpRLRL<SCALAR>(i, j, k, l, x) );
    _ops.push_back(temp_ptr);
  }

  // unknown keyword
  else
    throw pcps::Exception( (boost::format("error while reading hamiltonian input file.  \"%s\" is not a valid keyword") % key).str() );

}

///////////////////////////////////////////////////////////////////////////////////////////////////
/// \brief   loads the Hamiltonian from a text file
///
/// \param[in]   userinp  user input options
///
///////////////////////////////////////////////////////////////////////////////////////////////////

template <class SCALAR> void pcps::Hamiltonian<SCALAR>::initialize_from_file(const pcps::Input &userinp) {

  // get MPI info
  const MPI::Comm & comm = MPI::COMM_WORLD;
  const int nproc = comm.Get_size();
  const int myrank = comm.Get_rank();

  // allocate a buffer for storing input file lines
  const int buffer_size = 4096;
  boost::shared_array<char> buffer = pcps::allocate_shared_array<char>(buffer_size+1);

  // master process
  if (myrank == 0) {

    // open the input file
    const std::string infile_name = ( boost::format("%s%s") % userinp.workdir() % "pcps_hamiltonian.txt" ).str();
    std::ifstream input_file(infile_name.c_str());

    // check that the file opened successfully
    if (!input_file.is_open())
      throw pcps::Exception((boost::format("failed to open hamiltonian file \"%s\"") % infile_name).str());

    // get the first line of the input file
    input_file.getline(buffer.get(), buffer_size);

    // process the input file lines, which should contain a keyword followed by one or more values
    while (!input_file.eof()) {

      // create a string from the line and remove any comments (text after a %)
      std::string line(buffer.get());
      pcps::strip_comments(line, '%');

      // send the line to the slave processes
      int line_len = line.size();
      comm.Bcast((void *)(&line_len), 1, MPI::INT, 0);
      for (int j = 0; j < line_len; j++)
        buffer[j] = line.at(j);
      comm.Bcast((void *)buffer.get(), line_len, MPI::CHAR, 0);

      // add the operator described by the input line
      this->add_operator_from_file(line);

      // get the next input file line
      input_file.getline(buffer.get(), buffer_size);

    }

    // close the input file
    input_file.close();

    // signal the slave processes to stop receiving input data
    int stop_signal = -1;
    comm.Bcast((void *)(&stop_signal), 1, MPI::INT, 0);

  // slave processes
  } else {

    // receive the words in the input file lines read by the master process
    while (true) {

      // receive the length of the line
      int line_len;
      comm.Bcast((void *)(&line_len), 1, MPI::INT, 0);

      // if the stop signal was sent, exit the loop
      if (line_len < 0) break;

      // receive the line
      comm.Bcast((void *)buffer.get(), line_len, MPI::CHAR, 0);
      std::string line(buffer.get(), line_len);

      // add the operator described by the input line
      this->add_operator_from_file(line);

    }

  }

}

///////////////////////////////////////////////////////////////////////////////////////////////////
/// \brief   computes the modified one electron integrals
///
///   The quantum chemistry electronic hamiltonian can be written as
///
///   \f[ H = \sum_{i j} h_{i j} a^\dag_i a_j + \frac{1}{2} \sum_{ p q r s} [pq|rs] a^\dag_p a^\dag_r a_s a_q \f]
///
///   where \f$h_{i j}\f$ are the one electron integrals stored in _oei and
///   \f$[pq|rs]\f$ are the two electron integrals stored in _tei.
///
///   The local energy function in quantum monte carlo is more easily computed if we
///   rearrange the Hamiltonian as 
///
///   \f[ H = \sum_{i j} f_{i j} a^\dag_i a_j + \frac{1}{2} \sum_{ p q r s} [pq|rs] a^\dag_p a_q a^\dag_r a_s \f]
///
///   where we have commuted \f$a_q\f$ to be next to \f$a^\dag_p\f$ and absorbed the
///   result of the delta function with \f$a^\dag_r\f$ into the one electron integrals,
///   creating the modified one electron integrals
///
///   \f[  f_{i j} = h_{i j} - \frac{1}{2} \sum_p  [ip|pj]  \f]
///
///   These modified one electron integrals are stored in _mei.
///
/// \param[in]   userinp  user input options
///
///////////////////////////////////////////////////////////////////////////////////////////////////

template <class SCALAR> void pcps::Hamiltonian<SCALAR>::initialize_mei(const pcps::Input &userinp) {

  // modify the one electron integrals
  for (int i = 0; i < _norb; i++)
  for (int j = 0; j < _norb; j++) {
    _mei[i*_norb+j] = _oei[i*_norb+j];
    for (int p = 0; p < _norb; p++)
      _mei[i*_norb+j] -= 0.5 * this->get_tei(i,p,p,j);
  }

}

///////////////////////////////////////////////////////////////////////////////////////////////////
/// \brief   determines which pairs of orbitals have non-trivial two electron integrals
///
///   Finds each orbital index pair pq for which at least one of the integrals [pq|**]
///   is nonzero.  Adds all such pairs to the vector of nonzero pairs.
///
/// \param[in]   userinp  user input options
///
///////////////////////////////////////////////////////////////////////////////////////////////////

template <class SCALAR> void pcps::Hamiltonian<SCALAR>::initialize_nz_pairs(const pcps::Input &userinp) {

  // count the number of non-zero pairs
  int nnz = 0;
  for (int i = 0; i < _norb; i++)
  for (int j = 0; j < _norb; j++) {
    bool trivial = true;
    for (int k = 0; k < _norb; k++)
    for (int l = 0; l < _norb; l++)
      if ( std::abs(this->get_tei(i,j,k,l)) > userinp.tei_screening() )
        trivial = false;
    if (!trivial)
      nnz++;
  }

  // allocate space for the pairs
  _nz_pairs.assign(2*nnz, 0);

  // record the nonzero pairs
  std::vector<int>::iterator it = _nz_pairs.begin();
  for (int i = 0; i < _norb; i++)
  for (int j = 0; j < _norb; j++) {
    bool trivial = true;
    for (int k = 0; k < _norb; k++)
    for (int l = 0; l < _norb; l++)
      if ( std::abs(this->get_tei(i,j,k,l)) > userinp.tei_screening() )
        trivial = false;
    if (!trivial) {
      *it++ = i;
      *it++ = j;
    }
  }
  if ( it != _nz_pairs.end() )
    throw pcps::Exception("did not find the correct number of nonzero pairs pcps::Hamiltonian::initialize_nz_pairs");

  // get MPI info
  const MPI::Comm & comm = MPI::COMM_WORLD;
  const int nproc = comm.Get_size();
  const int myrank = comm.Get_rank();

  if (myrank == 0)
    std::cout << boost::format("found %i non-zero pairs in the two electron integrals") % nnz << std::endl << std::endl;

}

///////////////////////////////////////////////////////////////////////////////////////////////////
/// \brief   loads the Hamiltonian from PSI3
///
///   The full quantum chemistry hamiltonian is loaded using PSI3's molecular orbital integrals.
///   If local orbitals are desired, then PSI3's localize and transqt modules should be used to
///   convert the molecular orbitals to a local orbital basis.
///
/// \param[in]   userinp  user input options
///
///////////////////////////////////////////////////////////////////////////////////////////////////

template <class SCALAR> void pcps::Hamiltonian<SCALAR>::initialize_from_psi3(const pcps::Input &userinp) {

  // get MPI info
  const MPI::Comm & comm = MPI::COMM_WORLD;
  const int nproc = comm.Get_size();
  const int myrank = comm.Get_rank();

  // ensure that we are working with fermions
  if ( userinp.particle_type() != pcps::keyword::Part_Fermion )
    throw pcps::Exception( ( boost::format("particle_type must be \"%s\" in pcps::Hamiltonian::initialize_from_psi3")
                             % pcps::keyword::Part_Fermion ).str() );

  // allocate storage for the integrals
  _oei.assign(_norb*_norb, pcps::zero<SCALAR>());
  _mei.assign(_norb*_norb, pcps::zero<SCALAR>());
  _tei.assign( _norb*(_norb+1) * ( _norb*(_norb+1) + 2 ) / 8, pcps::zero<SCALAR>());

  // read integrals on root process
  if (myrank == 0)
    pcps::get_psi3_ints(userinp, _enuc, &_oei.at(0), &_tei.at(0));

  // send integrals to all processes
  comm.Barrier();
  comm.Bcast(&_enuc, 1, MPI::DOUBLE, 0);
  comm.Bcast(&_oei.at(0), _oei.size() * sizeof(SCALAR), MPI::CHAR, 0);
  comm.Bcast(&_tei.at(0), _tei.size() * sizeof(SCALAR), MPI::CHAR, 0);

//  // print nuclear repulsion energy
//  std::cout << boost::format("enuc = %20.12f") % _enuc << std::endl;
//  std::cout << std::endl;
//
//  // print one electron integrals
//  std::cout << boost::format("one electron integrals:") << std::endl;
//  std::cout << std::endl;
//  for (int i = 0; i < _norb; i++)
//  for (int j = 0; j < _norb; j++)
//    std::cout << boost::format(" %4i %4i %20.12f") % i % j % _oei[i*_norb+j] << std::endl;
//  std::cout << std::endl;
//
//  // print one electron integrals
//  std::cout << boost::format("two electron integrals:") << std::endl;
//  std::cout << std::endl;
//  for (int i = 0; i < _norb; i++)
//  for (int j = 0; j < _norb; j++)
//  for (int k = 0; k < _norb; k++)
//  for (int l = 0; l < _norb; l++)
//    std::cout << boost::format(" %4i %4i %4i %4i %20.12f") % i % j % k % l % this->get_tei(i,j,k,l) << std::endl;
//  std::cout << std::endl;

  // initialize the modified one electron integrals
  this->initialize_mei(userinp);

  // initialize list of non-zero pair indices
  this->initialize_nz_pairs(userinp);

  //throw pcps::Exception("Stopping after initializing Hamiltonian from PSI3");

}

///////////////////////////////////////////////////////////////////////////////////////////////////
/// \brief   loads the Hamiltonian from a Molpro FCI dump file
///
/// \param[in]   userinp  user input options
///
///////////////////////////////////////////////////////////////////////////////////////////////////

template <class SCALAR> void pcps::Hamiltonian<SCALAR>::initialize_from_molpro(const pcps::Input &userinp) {

  // get MPI info
  const MPI::Comm & comm = MPI::COMM_WORLD;
  const int nproc = comm.Get_size();
  const int myrank = comm.Get_rank();

  // ensure that we are working with fermions
  if ( userinp.particle_type() != pcps::keyword::Part_Fermion )
    throw pcps::Exception( ( boost::format("particle_type must be \"%s\" in pcps::Hamiltonian::initialize_from_molpro")
                             % pcps::keyword::Part_Fermion ).str() );

  // allocate storage for the integrals
  _oei.assign(_norb*_norb, pcps::zero<SCALAR>());
  _mei.assign(_norb*_norb, pcps::zero<SCALAR>());
  _tei.assign( _norb*(_norb+1) * ( _norb*(_norb+1) + 2 ) / 8, pcps::zero<SCALAR>());

  // read integrals on root process
  if (myrank == 0)
    pcps::get_molpro_ints(userinp, _enuc, &_oei.at(0), &_tei.at(0));

  // send integrals to all processes
  comm.Barrier();
  comm.Bcast(&_enuc, 1, MPI::DOUBLE, 0);
  comm.Bcast(&_oei.at(0), _oei.size() * sizeof(SCALAR), MPI::CHAR, 0);
  comm.Bcast(&_tei.at(0), _tei.size() * sizeof(SCALAR), MPI::CHAR, 0);

//  // print nuclear repulsion energy
//  std::cout << boost::format("enuc = %20.12f") % _enuc << std::endl;
//  std::cout << std::endl;
//
//  // print one electron integrals
//  std::cout << boost::format("one electron integrals:") << std::endl;
//  std::cout << std::endl;
//  for (int i = 0; i < _norb; i++)
//  for (int j = 0; j < _norb; j++)
//    std::cout << boost::format(" %4i %4i %20.12f") % i % j % _oei[i*_norb+j] << std::endl;
//  std::cout << std::endl;
//
//  // print one electron integrals
//  std::cout << boost::format("two electron integrals:") << std::endl;
//  std::cout << std::endl;
//  for (int i = 0; i < _norb; i++)
//  for (int j = 0; j < _norb; j++)
//  for (int k = 0; k < _norb; k++)
//  for (int l = 0; l < _norb; l++)
//    std::cout << boost::format(" %4i %4i %4i %4i %20.12f") % i % j % k % l % this->get_tei(i,j,k,l) << std::endl;
//  std::cout << std::endl;

  // initialize the modified one electron integrals
  this->initialize_mei(userinp);

  // initialize list of non-zero pair indices
  this->initialize_nz_pairs(userinp);

}

///////////////////////////////////////////////////////////////////////////////////////////////////
/// \brief   Loads the quantum chemistry Hamiltonian from text files.
///
/// \param[in]   userinp     user input options
/// \param[in]   order_1212  Indicates two electron integral order.  1212 if true, 1122 if false.
///
///////////////////////////////////////////////////////////////////////////////////////////////////

template <class SCALAR> void pcps::Hamiltonian<SCALAR>::initialize_qc_from_file(const pcps::Input &userinp, const bool order_1212) {

  // get MPI info
  const MPI::Comm & comm = MPI::COMM_WORLD;
  const int nproc = comm.Get_size();
  const int myrank = comm.Get_rank();

  // ensure that we are working with fermions
  if ( userinp.particle_type() != pcps::keyword::Part_Fermion )
    throw pcps::Exception("particle_type must be \"%s\" in pcps::Hamiltonian::initialize_qc_from_file") % pcps::keyword::Part_Fermion;

  // allocate storage for the integrals
  _oei.assign(_norb*_norb, pcps::zero<SCALAR>());
  _mei.assign(_norb*_norb, pcps::zero<SCALAR>());
  _tei.assign( _norb*(_norb+1) * ( _norb*(_norb+1) + 2 ) / 8, pcps::zero<SCALAR>());

  // read integrals on root process
  if (myrank == 0) {

    // read nuclear energy
    std::ifstream ne_ifs( (userinp.workdir() + "pcps_core_energy.txt").c_str() );
    if (!ne_ifs.is_open()) throw pcps::Exception("failed to open core energy file");
    ne_ifs >> _enuc;
    if ( ne_ifs.eof() )
      throw pcps::Exception("reached end of file while readng nuclear energy in pcps::Hamiltonian::initialize_qc_from_file");
    ne_ifs.close();

    // prepare some indices and a value for reading integrals
    int i, j, k, l;
    SCALAR val;

    // read one electron integrals
    std::ifstream oei_ifs( (userinp.workdir() + "pcps_oei.txt").c_str() );
    if (!oei_ifs.is_open()) throw pcps::Exception("failed to open one electron integral file");
    oei_ifs >> i >> j >> val;
    while (!oei_ifs.eof()) {
      if ( i < 0 || i >= _norb ) throw pcps::Exception("bad index when reading one electron integrals: %i") % i;
      if ( j < 0 || j >= _norb ) throw pcps::Exception("bad index when reading one electron integrals: %i") % j;
      _oei[i*_norb+j] = _oei[j*_norb+i] = val;
      oei_ifs >> i >> j >> val;
    }
    oei_ifs.close();

    // read two electron integrals
    std::ifstream tei_ifs( (userinp.workdir() + "pcps_tei.txt").c_str() );
    if (!tei_ifs.is_open()) throw pcps::Exception("failed to open two electron integral file");
    tei_ifs >> i >> j >> k >> l >> val;
    while (!tei_ifs.eof()) {
      if ( i < 0 || i >= _norb ) throw pcps::Exception("bad index when reading two electron integrals: %i") % i;
      if ( j < 0 || j >= _norb ) throw pcps::Exception("bad index when reading two electron integrals: %i") % j;
      if ( k < 0 || k >= _norb ) throw pcps::Exception("bad index when reading two electron integrals: %i") % k;
      if ( l < 0 || l >= _norb ) throw pcps::Exception("bad index when reading two electron integrals: %i") % l;
      ( order_1212 ? this->get_tei(i,k,j,l) : this->get_tei(i,j,k,l) ) = val;
      tei_ifs >> i >> j >> k >> l >> val;
    }
    tei_ifs.close();

  }

  // send integrals to all processes
  comm.Barrier();
  comm.Bcast(&_enuc, 1, MPI::DOUBLE, 0);
  comm.Bcast(&_oei.at(0), _oei.size() * sizeof(SCALAR), MPI::CHAR, 0);
  comm.Bcast(&_tei.at(0), _tei.size() * sizeof(SCALAR), MPI::CHAR, 0);

  // initialize the modified one electron integrals
  this->initialize_mei(userinp);

  // initialize list of non-zero pair indices
  this->initialize_nz_pairs(userinp);

}

///////////////////////////////////////////////////////////////////////////////////////////////////
/// \brief   constructs the Hamiltonian according to the user's input options
///
/// \param[in]   userinp  user input options
///
///////////////////////////////////////////////////////////////////////////////////////////////////

template<class SCALAR> pcps::Hamiltonian<SCALAR>::Hamiltonian(const pcps::Input &userinp)
  : _norb(userinp.nsites()/2),
    _enuc(0.0),
    _qc(    userinp.hamiltonian_type() == pcps::keyword::Ham_PSI3
         || userinp.hamiltonian_type() == pcps::keyword::Ham_Molpro
         || userinp.hamiltonian_type() == pcps::keyword::Ham_Text_1212 
         || userinp.hamiltonian_type() == pcps::keyword::Ham_Text_1122 )
{

  // get MPI info
  const MPI::Comm & comm = MPI::COMM_WORLD;
  const int nproc = comm.Get_size();
  const int myrank = comm.Get_rank();

  // setup for the Hubbard Hamiltonian
  if ( userinp.hamiltonian_type() == pcps::keyword::Ham_Hubbard )
    this->initialize_hubbard(userinp);

  // setup for the Heisenberg Hamiltonian
  else if ( userinp.hamiltonian_type() == pcps::keyword::Ham_Heisenberg )
    this->initialize_heisenberg(userinp);

  // setup for a user inputted Hamiltonian
  else if ( userinp.hamiltonian_type() == pcps::keyword::Ham_File )
    this->initialize_from_file(userinp);

  // setup for a quantum chemistry hamiltonian from PSI3
  else if ( userinp.hamiltonian_type() == pcps::keyword::Ham_PSI3 )
    this->initialize_from_psi3(userinp);

  // setup for a quantum chemistry hamiltonian from Molpro
  else if ( userinp.hamiltonian_type() == pcps::keyword::Ham_Molpro )
    this->initialize_from_molpro(userinp);

  // setup for a quantum chemistry hamiltonian from text files in 1212 ordering
  else if ( userinp.hamiltonian_type() == pcps::keyword::Ham_Text_1212 )
    this->initialize_qc_from_file(userinp, true);

  // setup for a quantum chemistry hamiltonian from text files in 1122 ordering
  else if ( userinp.hamiltonian_type() == pcps::keyword::Ham_Text_1122 )
    this->initialize_qc_from_file(userinp, false);

  // throw an error if the Hamiltonian type is unknown
  else 
    throw pcps::Exception( (boost::format("\"%s\" is not a valid hamiltonian type") % userinp.hamiltonian_type()).str() );

  // separate operators into diagonal and off-diagonal
  for (int i = 0; i < _ops.size(); i++)
    if ( typeid(*_ops.at(i)) == typeid(pcps::OpRL<SCALAR>) || typeid(*_ops.at(i)) == typeid(pcps::OpRLRL<SCALAR>) )
      _ops_off_diag.push_back(_ops.at(i));
    else
      _ops_diag.push_back(_ops.at(i));

  // separate operators by type
  for (typename std::vector<typename pcps::OpBase<SCALAR>::const_ptr>::const_iterator op = _ops.begin(); op != _ops.end(); op++) {
    if ( typeid(**op) == typeid(pcps::OpRL<SCALAR>) ) {
      boost::shared_ptr<const pcps::OpRL<SCALAR> > derived_op_ptr
        = boost::dynamic_pointer_cast<const pcps::OpRL<SCALAR>, const pcps::OpBase<SCALAR> >(*op);
      _opsRL.push_back(*derived_op_ptr);
    } else if ( typeid(**op) == typeid(pcps::OpRLRL<SCALAR>) ) {
      boost::shared_ptr<const pcps::OpRLRL<SCALAR> > derived_op_ptr
        = boost::dynamic_pointer_cast<const pcps::OpRLRL<SCALAR>, const pcps::OpBase<SCALAR> >(*op);
      _opsRLRL.push_back(*derived_op_ptr);
    } else if ( typeid(**op) == typeid(pcps::OpNN<SCALAR>) ) {
      boost::shared_ptr<const pcps::OpNN<SCALAR> > derived_op_ptr
        = boost::dynamic_pointer_cast<const pcps::OpNN<SCALAR>, const pcps::OpBase<SCALAR> >(*op);
      _opsNN.push_back(*derived_op_ptr);
    } else if ( typeid(**op) == typeid(pcps::OpN<SCALAR>) ) {
      boost::shared_ptr<const pcps::OpN<SCALAR> > derived_op_ptr
        = boost::dynamic_pointer_cast<const pcps::OpN<SCALAR>, const pcps::OpBase<SCALAR> >(*op);
      _opsN.push_back(*derived_op_ptr);
    } else if ( typeid(**op) == typeid(pcps::OpIdentity<SCALAR>) ) {
      boost::shared_ptr<const pcps::OpIdentity<SCALAR> > derived_op_ptr
        = boost::dynamic_pointer_cast<const pcps::OpIdentity<SCALAR>, const pcps::OpBase<SCALAR> >(*op);
      _opsIdentity.push_back(*derived_op_ptr);
    } else
      throw pcps::Exception( (boost::format("unexpected operator type in Hamiltonian constructor: \"%s\"") % typeid(**op).name()).str() );
  }

  // print some information about the hamiltonian
  if ( myrank == 0 && !_qc ) {

    std::cout << "hamiltonian contains " << _ops.size() << " operators" << std::endl;
    std::cout << std::endl;

    for (typename std::vector<OpBaseConstPtr>::const_iterator op = _ops.begin(); op != _ops.end(); op++) {
      std::cout << boost::format("%40s") % pcps::type_lookup(typeid(**op));
      std::cout << boost::format("  %s") % pcps::format_number("%10.4f", (**op).coeff());
      for (int i = 0; i < (**op).sites().size(); i++)
        std::cout << boost::format(" %4i") % (**op).sites().at(i);
      std::cout << std::endl;
    }
    std::cout << std::endl;

  } else if ( myrank == 0 && _qc ) {

    std::cout << boost::format("nuclear and core energy = %20.12f") % _enuc << std::endl;
    std::cout << std::endl;

  }

}

// explicitly instantiate the constructor
template pcps::Hamiltonian< double               >::Hamiltonian(const pcps::Input &);
template pcps::Hamiltonian< std::complex<double> >::Hamiltonian(const pcps::Input &);
